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Abstract 

We present a diffusion Monte Carlo study of a single vortex in two- 
dimensional superfluid liquid 4 He within the fixed node approximation. We 
use both the Feynman phase and an improved phase which includes backflow 
correlations to model the nodal surface of the vortex wavefunction. Results 

for the particle density, core radius and excitation energies are presented. 
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Diffusion and Green's function Monte Carlo simulations have become a standard tool in 
the study of ground-state properties of Bose quantum liquids at zero temperature. However, 
the development of ab initio computational methods to investigate properties of excited 
states, such as the phonon-roton branch and vortex excitations in superfluid 4 He, is still 
a challenging problem at the forefront of present research in the field of computational 
techniques applied to condensed matter physics. In this Letter we present results on the 
microscopic structure of a single vortex excitation in two-dimensional liquid 4 He at zero 
temperature, obtained by employing a fixed-node diffusion Monte Carlo (DMC) method [|TJ. 

The idea that circulation in superfluid 4 He is quantized is due to Onsager 0, whereas 
the possibility that vortices might have a core of atomic dimensions was first put forward 
by Feynman || who also proposed a microscopic wavefunction for the vortex state. Gross, 
Pitaevskii and Fetter [|J investigated the structure of vortex states in a weakly interacting 
inhomogeneous Bose gas using a mean field approach. The first attempt to study quan- 
tized vortices in a strongly interacting system, such as liquid 4 He, is due to Chester, Metz 
and Reatto || who calculated the energy of a vortex line with the use of a variational ap- 
proach involving integral equations. Only quite recently have appeared new calculations 
of the structure of vortex states in superfluid 4 He |3|-[|. The very recent paper by Ortiz 
and Ceperley |J is the first attempt to tackle the problem of the vortex core structure by 
employing ab initio computational techniques. Our method in the present Letter is similar 
to that used by these authors, but the approach is different and our results for the core 
energy and the particle density near the vortex axis are significantly different from the ones 
obtained in Ref. |J. 

A vortex excitation is an eigenstate of the iV-particle Hamiltonian H and of the z- 
component of the angular momentum L z with eigenvalue hN£, corresponding to an integer 
number I of quanta of circulation [|Kj . The simplest microscopic wavefunction to describe 
a vortex state was introduced by Feynman [[J: 
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where ifF = Z)t=i n @i is the Feynman phase with d{ the azimuthal angle of the z-th particle, 
$q(R) describes the ground state of the system and /(r$) is a function of the radial distance 
of each particle from the vortex axis, which models the density near the core. In what follows 
we only consider vortices with one quantum of circulation, i.e., i = ±1. In the recent paper 
by Ortiz and Ceperley |J, a systematic method to improve the phase of the wavefunction 
is devised. Starting from the Feynman phase ^ as zeroth order ansatz, the first correction 
includes backflow correlations giving an improved phase of the form 

(pBF = V>F + ><22'y{ri,rj,rij)—^ii($i-9 j ) . (2) 

i¥=3 r% 

The wavefunction constructed with the phase <£bf is the vortex analogue of the Feynman- 



Cohen backflow wavefunction for the phonon-roton excitation branch fLl |. 

To go beyond a variational estimation of the properties of the vortex state, given by the 
above model wavef unctions, we have used a DMC method. This method solves the many- 
body Schrodinger equation in imaginary time for the function /(R,i) = ?/V(R)$(R, t) 

df(K,t) 



dt 



-DVy(R, t) + DVr(F(R)/(R, t)) + (E L (R) - E)f(R, t) , (3) 



where $(R, t) is the wavefunction of the system and ^t(R) is a trial function used for 
importance sampling. In the above equation, El(H) = ift^ (R)-^V't(R) is the local energy 
and F(R) = 2^ 7 ; 1 (R)Vr^t(R) is the so-called quantum drift force; D = h 2 /2m, with m the 
mass of the particles, plays the role of a diffusion constant, R stands for the 3iV-coordinate 
vector of the iV particles of the system and E is an arbitrary energy shift. Equation ([3D 
is a diffusion equation for the probability distribution /(R, t) which evolves in time due to 
diffusion, drift and branching processes. If $(R) represents the wavefunction of the lowest 
energy eigenstate of the system not orthogonal to the trial function ip?, the asymptotic 
solution of Eq. (^) is given by /(R, t — ► oo) = ^r(R)$(R) and the corresponding energy 
eigenvalue can be calculated exactly. 

In order to deal with a real walker probability distribution function /(R, t), we choose as 
trial function the superposition of two vortex states, one with positive and one with negative 



circulation, which are degenerate in energy. The importance sampling function has therefore 
the form 

N 

^ T (R) = cos(^(R))n/(^r(R) > (4) 

i=i 

where ip^ is the importance sampling function for the ground state and for the phase p 
we have used both the Feynman phase p>F and the backflow phase <pbf- The sign problem 
associated with the use of the above trial wavefunction in the DMC algorithm has been 
dealt with in the framework of the fixed- node (FN) approximation Q. This approach, 
which has been extensively used in the calculation of ground-state properties of fermionic 
systems, yields an upper bound to the energy eigenvalue |]. The quantum drift force acting 
on each particle, as obtained from the trial wavefunction (|j), can be written as the sum 
F*(R) = F^(R) + F?>(R). The first term in the sum is independent of the phase ip of the 
wavefunction, whereas the second term contains the ip dependence and has the form 

F*(R) = -2tan(^)Vi^ , (5) 

where p>i is the contribution of the z-th particle to the collective phase, <p = J2i=i,N <Pi- m 
the same way, the local energy £x(R) can be decomposed in the sum of a phase independent 
term En and a phase dependent term El 2 which is given by 

N t 1 \ 

E L2 (R) =DY, ((V^) 2 + tan(^)V 2 ^ - -F\(R) ■ F*(R)J . (6) 

Vortex states are characteristic of systems with rotational invariance, but, at the same 
time, simulations must deal with a finite number of particles. A simple choice is to restrict 
the N particles to be inside a cylindric box with the vortex at the center and rigid boundary 
conditions on the walls. This is the geometry chosen for example in Ref. ||. By confining 
the particles some problems arise, such as the choice of the confining potential and surface 
effects which can be relevant if the box size is not large enough. In the present work, we have 
removed the confinement in order to keep surface effects as small as possible. An important 
point one needs to solve when attempting to simulate a vortex excitation in an infinite 
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system, is the calculation of the collective phase <^(R). For a very large system one expects 
that the features of the vortex state near the core are weakly influenced by the behaviour 
of particles which are far from the vortex axis. It is thus reasonable to assume that, if the 
collective phase is decomposed in the sum 

where the first term sums the contributions to the phase of all the particles with distance 
from the vortex axis within the cutoff length f and the second term gives the contribution 
coming from all the other particles, the phase fluctuations in the second term are irrelevant 
for the core structure of the vortex and it can be safely approximated by its mean value 
J^nyrft — 0. If this prescription is employed in the calculation of the collective phase 
entering the expressions (|5|) and (|6]) for the drift force and the local energy, one expects 
that for a large enough cutoff length f the properties of the vortex state near the axis are 
properly simulated. However, the straightforward interpretation of the decomposition (|7|) 
would be of no practical use because the collective phase would change discontinously by 
a large amount each time a particle exits or enters the region delimited by the cutoff f. 
Instead, the procedure we have adopted is to use the decomposition (|7|) as a way of tagging 
the particles that will contribute to the collective phase for a long simulation run. The 
dependence of the results on the cutoff length have been studied and as will be discussed 
later no appreciable changes are seen for values of f larger than approximately 3 a (a = 2.556 
A). It is worth noticing that in the limit f —> the collective phase vanishes and all the 
terms containing explicitly ip in the expressions fl5|) and (^j for the drift force and the local 
energy disappear. In this case, the DMC calculation can be shown to be equivalent to 
the fixed-phase (FP) method employed by Ortiz and Ceperley in Ref. 0, where the phase 
of the wavefunction is fixed and the DMC algorithm is used to solve the equation for the 
modulus of the wavefunction. By taking a finite value for the cutoff f, one allows for phase 
fluctuations in the system around the phase introduced with the trial function, and in the 
limit f — > oo one recovers the full FN approximation. Our results actually show that a finite 
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cutoff is enough to account completely for the phase fluctuations. 

Once established that the contribution of the more distant particles to the collective phase 
is irrelevant, there is no a compelling reason for not extending the system using periodic 
boundary conditions. Surface effects exist also in this case, in the sense that particles near 
the walls of the box "see" the artificial density perturbations associated with the image 
vortices in the adjacent boxes. However, these effects are negligible for a reasonable size of 
the simulation cell and are definitely smaller than the density oscillations induced by rigid 
walls. 

We are now in a position to discuss our results. We consider 64 particles in a two- 
dimensional box of length L = 15.0 a, which corresponds approximately to the equilibrium 
density of the two-dimensional homogeneous liquid ||12|| . The atoms interact through the two- 
body HFD-B(HE) potential which is a revised version of the HFDHE2 Aziz potential. 
The ground-state trial wavefunction we have chosen is the McMillan two body function 

(R) = ]lj<i ex P(~ b 5 /2rfj) with b = 1.205 o as in the ground-state calculation of Ref. []12 



For the radial function /(r), which models the structure of the vortex core, we consider two 
different options: /i(r) = 1 — e^~ r ^ a \ and /2( r ) = 1- The first function gives a density in 
the trial function which decreases to zero at the vortex axis over a distance of order a, for 
which we take the value a — 1 A, whereas the second one does not contain any parameter 
associated with the vortex core. For the backflow function 7 entering Eq. (JJ) we have used 
the same functional form 7(7-3, rj, r^) = exp(— a(rf + rj) — f3rf-) and the same values for 
the parameters a, j3 and A as in Ref. ||. We have performed the calculation using three 
different trial wavefunctions: ipr 1 an d V'r 2 which correspond to the Feynman phase tpp 
with the radial terms fi and f 2 respectively, and ip^ corresponding to the backflow phase 
ifBF with f 2 . In Fig. 1 we show results for the particle density p(r) using ip^ 1 , ip^ 2 and 
ipT ■ These results have been obtained by means of mixed estimators which are significantly 
biased by the choice of the trial wavefunction. As apparent from Fig. 1, the behaviour near 
the core is strongly affected by the introduction or not of the radial term fi{r). In order to 
remove the influence of the trial wavefunction in mixed estimators of coordinate operators 



one can use pure estimators. In the present work we have employed the method presented in 



Ref . [14j] . The results for the pure density profiles do not depend anymore on the use or not 
of the radial term fi(r) and the three cases converge to a single result. The common pure 
profile is presented in Fig. 2, where it clearly appears that a zero particle density is reached 
on the vortex axis. This result is in contrast with the prediction of a significant non-zero 
particle density on the axis obtained with the backflow phase in Ref. ||. In our opinion the 
result of Ref. || can be influenced by the extrapolation technique used by the authors to 
improve the mixed estimator result. In fact, the usual linear extrapolation technique, which 
is accurate to the same order as the one employed in Ref. |J, would change their result for 
the particle density on the vortex axis in an amount comparable to their prediction. 

In Table I are reported the energies per particle for the different trial wavef unctions, 
together with the energy per particle E /N in the ground state. The two results for the 
Feynman phase are almost equal, whereas in the case of the backflow phase the system is 
slightly more bounded in accordance with the improvement of the nodal surface. The cutoff 
length f for the calculation of the collective phase has been taken as f = 6 a. In the case 
of the Feynman phase reducing the cutoff length does not give any change in the results 
down to the FP limit f = 0. For the backflow phase, though, the FP energy is slightly less 
negative E^/N = -0.8136 ± 0.0020 K. 

The excitation energy E v {r) of a vortex inside a disk of radius r is obtained as the differ- 
ence between the total energy of the disk with and without the vortex E v (r) = E(r) —E (r). 
For large distances from the vortex axis E v (r) is usually decomposed in an hydrodynamic tail, 
which depends logarithmically on r, and a core energy E c : E v {r) = nh 2 po/m log(r/£) + E c , 
where po is the homogeneous density far from the axis and £ is the vortex core radius. In Fig. 
3 we show the vortex excitation energy E v (r) for the different choices of the trial wavefunc- 
tion. For distances r > 6 A E v (r) shows the expected hydrodynamic behaviour with a small 
negative shift of the backflow energy with respect to the Feynman one in agreement with 
the total energy results shown in Table I. For small r's, the estimation of E v (r) is not exact 
and exhibits the influence of the trial wavefunction used. It is worth noticing the absence of 



spurious oscillations at large r, present in the results of Ref. 0, which are induced by the 
confining rigid walls. 

The core radius £ can be estimated, following Ref. ||, as the position of the maximum in 
the azimuthal circulating current J$(r). The radial dependence of Je(r) has been estimated 
from the pure density profile using the expression for the current in the FP approximation 
at the Feynman level, Jq(t) = p p (r)/r. The value obtained is £ = 2.10 ± 0.20 A which is in 
agreement with the result reported in Ref. |J. By calculating the hydrodynamic contribution 
to the total energy for our square simulation box, we get the core energy 1 = 1.23 ± 0.25 
K, EP = 1.18 ± 0.26 K and E^ F = 1.00 ± 0.26 K for the Feynman and backflow phases, 
respectively. These values coincide with the results obtained by a fit to E v (r) for r > 6 A. 
Our results for E c are significantly smaller from the ones obtained in Ref. || and the values 
of E F are close to the variational results of Ref. || based on the Feynman phase. 

In conclusion, we have studied the structure of a vortex excitation in two-dimensional 
superfluid 4 He using a DMC method within the fixed node approximation. The collective 
phase of the vortex has been dealt with a method that allows for the use of periodic boundary 
conditions, removing spurious surface effects introduced by the use of confining geometries. 
The fixed phase approximation is recovered as a limiting case in our approach. The result 
for the density profile predicts a zero particle density on the vortex axis for the two model 
phases used. On the other hand, the inclusion of backflow correlations in the phase gives a 
slightly smaller upper bound for the excitation energy. Finally, we would mention that the 
fixed node DMC method used in the present work permits the study of other excited states 
in liquid 4 He. Work is in progress to extend our calculations to the phonon-roton branch 
and the vortex-antivortex excitation in two dimensions. 
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FIGURES 

FIG. 1. Mixed density profiles. Solid line, with ip F1 ; short-dashed line, with if) F2 , and 
long-dashed line with tp^ F . 



FIG. 2. Pure density profile. 

FIG. 3. Radial dependence of the excitation energy. Solid line, with if)? 1 ', short-dashed line, 
with i/j F2 , and long-dashed line with ip^ F . 
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TABLES 

TABLE I. Results for the energies per particle. Epi/N, Ep2/N, and Ebf/N correspond to 
the trial functions V't 1 ' V't 2 an d V't^; respectively. Eq/N is the ground-state energy. 



E Fl /N(K) 


E F2 /N(K) 


E BF /N(K) 


Eq/N(K) 


-0.8162(16) 


-0.8171(18) 


-0.8199(18) 


-0.8957(25) 
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